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We derive expressions of interatomic force and heat current for many-body potentials such as 
the Tersoff, the Brenner, and the Stillinger-Weber potential used extensively in molecular dynamics 
simulations of covalently bonded materials. Although these potentials have a many-body nature, 
a pairwise force expression that follows Newton’s third law can be found without referring to any 
partition of the potential. Based on this force formula, a stress applicable for periodic systems can 
be unambiguously defined. The force formula can then be used to derive the heat current formulas 
using a natural potential partitioning. Our heat current formulation is found to be equivalent to 
most of the seemingly different heat current formulas used in the literature, but to deviate from the 
stress-based formula derived from two-body potential. We validate our formulation numerically on 
various systems descried by the Tersoff potential, namely three-dimensional silicon and diamond, 
two-dimensional graphene, and quasi-one-dimensional carbon nanotube. The effects of cell size and 
time used in the simulation are examined. 

PACS numbers: 02.70.Ns, 05.60.Cd, 44.10.-|-i, 66.70.-f 

I. INTRODUCTION 

Molecular dynamics (MD) simulation has been used 
extensively to study thermal transport properties of ma¬ 
terials. There are mainly two methods for computing lat¬ 
tice thermal conductivity in the level of classical MD sim¬ 
ulations: the direct method^i^ [also called the nonequi¬ 
librium MD (NEMD) method] based on the Fourier’s law 
and the Green-Kubo^4^ method (also called the equilib¬ 
rium MD method) based on the Green-Kubo formula. 

Gross-checking of these two methods has also been the 
subject of several worksi^— In the direct method, the 
thermal conductivity is usually computed by measuring 
the steady-state temperature gradient at a fixed external 
heat current, analogous to the experimental situation. 

In contrast, in the Green-Kubo method, the thermal con¬ 
ductivity is computed by integrating the heat current au¬ 
tocorrelation function (HCACF) using the Green-Kubo 
formula. While the heat current in the direct method is 
created by scaling the velocities in the source and sink 
regions of the simulated system, which does not depend 
on the underline interatomic potential, the heat current 
in the Green-Kubo method is the summation of the mi¬ 
croscopic heat currents of the individual atoms in the 
simulated system, which generally depends on the spe¬ 
cific interatomic potential used. 

For a two-body potential, where a pairwise force can be 
directly defined, the heat current expression used in the 


Green-Kubo formula is well established. It is currently 
implemented in Large-scale Atomic/Molecular Massively 
Parallel Simulator (LAMMPS)^ in terms of the per-atom 
stress and works well for systems described by two-body 
potentials such as Lennard-Jones argon. However, it is 
not widely recognized that the heat current expression 
based on the per-atom stress is only applicable to two- 
body potentials, and is not guaranteed to produce cor¬ 
rect results for systems described by a many-body poten¬ 
tial, such as the widely used Tersoff potentialjiS, Brenner 
potentialfii and Stillinger-Weber potentialJ^ In the lit¬ 
erature, there have been quite a few formulationsH^il of 
the heat current for the Tersoff/Brenner potential, which 
seem to be inequivalent to each other 

In this work, we present detailed derivations of the 
heat current expressions for these many-body potentials. 
We show that many of the seemingly different formula¬ 
tions of the heat current are equivalent, except for some 
marginal differences resulting from a different decomposi¬ 
tion of the total potential into site (per-atom) potentials. 
Our derivation is facilitated by establishing the existence 
of a pairwise force respecting Newton’s third law, which is 
not widely recognized so far. Based on the pairwise force, 
a well defined expression for the virial tensor can also be 
obtained. By comparing with finite-difference calcula¬ 
tions, we validate the proposed pairwise force expression 
unambiguously. The derived expression is equivalent to 
other alternatives which do not respect Newtons third 
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law explicitly, but it has an advantage of allowing for 
an efficient implementation on graphics processing units 
(GPUs), which attains a speedup factor of two orders of 
magnitude (compared to our optimized serial CPU code) 
for large simulation cell sizes. 

Using the efficient GPU code, we perform a com¬ 
prehensive validation of our formulations by calculating 
lattice thermal conductivities of various kinds of mate¬ 
rial described by the Tersoff potential, including three- 
dimensional (3D) silicon and diamond, two-dimensional 
(2D) graphene, and quasi-one-dimensional (QID) carbon 
nanotube (CNT). For each material, we examine the con¬ 
vergence of the calculated thermal conductivity with re¬ 
spect to the total simulation time, the correlation time, 
and the finite-size effects, before comparing our results 
with previous ones. Last, we present explicit numerical 
evidence that the stress-based heat current expression is 
inequivalent to our formulation for the Tersoff potential. 


II. THEORY 

A. Green-Kubo method for thermal conductivity 
calculations 

The Green-Kubo formula for the running thermal con¬ 
ductivity (RTC) tensor K^^(t) {fJL^v = x,y,z) at a given 
correlation time t can be expressed as^“— 

ftMt'(t) = (1) 

where ks is Boltzmann’s constant, T is the absolute tem¬ 
perature, and V is the volume of the simulation cell. The 
HCACF C^v{t) is defined as 

= = ( 2 ) 

the ensemble average of the product of two heat currents 
separated by t. In the MD simulation, the ensemble av¬ 
erage is substituted by a time average. The simulation 
time required for achieving high statistical accuracy of 
the computed thermal conductivity in the Green-Kubo 
method is usually quite challenging, as we show later. 
The Green-Kubo method is capable of calculating the 
full conductivity tensor, but the following cases are suffi¬ 
cient to verify our formulations: (1) isotropic 3D systems, 
such as diamond, where we define the conductivity scalar 
as {kxx + Kyy + Kzz)l‘i, (2) Isotroplc 2D systems, such as 
graphene, where we define the in-plane conductivity as 
(Ka;a;-|-/Cyy)/2, and (3) QID systems, such as CNT, where 
only the conductivity along the tube is needed. Periodic 
boundary conditions are needed in all the transport di¬ 
rections. In the following, we use J to represent the heat 
current vector with components J^, Jy, and Jz- 


B. General expression of the heat current 

The heat current used in Eq. ([2]) is defined as the time 
derivative of the sum of the moments of the site energies 

Ei = + Ui (3) 

of the particles in the system^: 

^ IE = E +E ^4^- 

i i i 

Here m^, v^, and Ui are the mass, velocity, and potential 
energy of particle i, respectively. Conventionally, one 
defines a kinetic part 

kin ~ ^ ^ '^iEi (5) 

i 

and a potential part 

i 

and write the total heat current as a sum of them: 

J — J kin “f J pot ■ (7) 

The kinetic term J^in needs no further derivation, apart 
from a possible issue of defining Ui for a many-body po¬ 
tential, and the potential term Jpot can be written as 

Jvot = ^r,{F,-Vi) + ^r~, ( 8 ) 

i i 

where the kinetic energy theorem, {^rmvf) = Fi ■ Vi, 
Fi being the total force on particle i, has been used. 
The kinetic term is also called the convective term, and 
is mostly important for gases. For Lennard-Jones liquid, 
Vogelsang et al^ showed that the thermal conductivity 
is mainly contributed by the partial HCACF involving 
the potential-potential term. For solids, the kinetic term 
barely contributes and can be simply discarded. Note 
that the kinetic and potential terms defined here cor¬ 
respond to the potential and kinetic terms, respectively, 
used in the Einstein formalism studied by Kinaci et al.^ 
who also found that the convective term (the potential 
term in the Einstein formalism) does not contribute to 
the thermal conductivity for solids. We thus focus on the 
potential part [Eq. ([5])] in the following discussions. 

C. Heat current for two-body potentials 

Before discussing many-body potentials, let us first ex¬ 
amine the case of two-body potentials. Eor these, the 
total potential energy of the system can be written as 

^ = ^EEt^b, (9) 

i 
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where the pair potential between particles i and j, Uij = 
Uji = Uij (vij ), only depends on the distance between 
the particles. The factor of 1/2 in the above equation 
compensates the double-counting of the pair potentials; 
one can equally omit it by requiring j > i (or j < i). The 
derived forces are purely pairwise and Newton’s third law 
is apparently valid: 


The potential part of the heat current is also intimately 
related to the virial part of the stress tensor. To see 
this, we first note that the virial W can be written as a 
summation of individual terms, 

W = ^W„ (18) 


F, (10) 

3¥=i 


^ dUij_ ^ _ 
dr,, ^ 


( 11 ) 


where Fij is the force on particle i due to particle j and 
the convention;^ 


( 12 ) 


where the per-atom virial for a periodic system reads 
'W, = -]-^r,j®F,j. (19) 

Therefore, the potential part of the heat current can be 
expressed in terms of the per-atom virial as 

jpor = (20) 

i 


for the relative position between two particles is adopted. 
If periodic boundary conditions are applied in a given 
direction, the minimum image convention is used to all 
the relative positions in that direction. Using the above 
notations, the first term on the right hand side (RHS) of 
Eq. (I5|) can be written as 

'^r,{F,-Vi)='^'^r,[F,j-Vi). (13) 

i i jiii 

To make further derivation for the second term on the 
RHS of Eq. ([5]), one has to make a choice for the site 
potential Ui. A natural choice is Ui = ^ ^ 

for two-body potentials, it does not matter much how to 
define the site potential. For example, the above choice is 
equivalent io Ui = \ '^j^iiUij + Uji) because Uij = Uji. 
Therefore, the second term on the RHS of Eq. ([S]) can be 
written as 

= (14) 

i i j^i 

Using the above two expressions, we can write the poten¬ 
tial term of the heat current as 

EE r,[F,j ■ {Vi + Vj)]. (15) 

i j^i 

In numerical calculations, the absolute positions, r^, will 
cause problems for systems with periodic boundary con¬ 
ditions. Fortunately, one can circumvent the difficulty 
by using the Newton’s third law Eq. (HU), from which we 
have 

Jlot = -^EE^bl^b ■ + (16) 

i j^i 

where only the relative positions, r^, are involved. This 
expression is also equivalent to a less symmetric form: 

= ( 17 ) 

i j^i 


The current implementation of the Green-Kubo formula 
for thermal conductivity in LAMMPS adopts this stress- 
based formula. However, as we show later, it does not 
apply to many-body potentials. 

D. Force expressions for Tersoff potential 

We now move on to many-body potentials, first focus¬ 
ing on the Tersoff potential. The total potential energy 
for a system described by the Tersoff potential can also be 
written in the form of Eq. with the “pair potential” 
Uij taking the form ofi^ 


Uij = fc{rij) [fR{r,j) - bijfA{rij)] , 

(21) 


(22) 

Cij — y ] fc{rik)gijk, 

(23) 

(P + {h — cos Oijk)"^ ' 

(24) 

Here, (3, n, c, d, and h are parameters and dyfc 
angle formed by and r^fc, which means that 

is the 

COS — COS Oif^j — 

rijVik 

(25) 

For simplicity, the dependence of the parameters on the 
particle type is omitted in the above equations. Detailed 
rules for determining the parameters in systems with two 
kinds of atom can be found in Ref. [^. While the 
functions fc, Jr, and Ja only depend on r^ , the bond- 
order function bij also depends on the positions of the 


neighbor particles of i and j and thus generally, Uij ^ 
Uji, which is a manifestation of the many-body nature of 
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the Tersoff potential. However, we notice that bij, hence 
Uij, is only a function of the position difference vectors 
originating from particle i (In the equation below, k = j 
is allowed.): 


Uij = Uij {{rik}k^i). (26) 


This property will play a crucial role in the following 
derivations. 

We now start to derive the force expressions for the 
Tersoff potential. We begin with the definition: 


F- = 

^ I - 


dU 

dri 


1 dUjk 

2 ^ ^ dri ' 

J k=^J 


(27) 


We can expand it as 


and Newton’s third law 

F Tersoff 171 Tersoff /nA\ 

ij ^ ji 

still holds. 

In the above derivations, we have not assumed any 
form of the site potential Ui. The definition of Ui for a 
many-body potential amounts to a decomposition of the 
total potential into site potentials. While such a decom¬ 
position is not needed for the derivation of the forces, it 
is needed for deriving the heat current, which involves 
time-derivative of the site potential [cf. Eq. ([5])]. A nat¬ 
ural choice for the decomposition is 

U = '^Ui with ^ X] (35) 

i j^i 


F 

I 


1 

2 


E dUik dUji 

dri ^ dri 

k^ti j^i 


+ EE 


dUjk 

dri 


(28) 

The first, second, and third terms on the RHS of Eq. (1251) 
correspond to the parts with j = i, k = i, and j^k ^ i in 
Eq. (|27ll . respectively. Then, using Eq. (l26l) . we have 


There is no clear physical intuition favoring this decom¬ 
position over others [cf. Eq. (IB18I) ]. but we find that 
Eq. (1351) is a very reasonable definition. To show this, we 
notice that the site potential defined by Eq. (1551) is also 
only a function of the relative positions originating from 
particle v. 


= ({<•«),V.). (36) 


EE 


dUik dr ij dUji dr jk 

... ... drij dri ^.^.dr^k dri 


1 dUjk drjm 

O Z -J / -J ■ 


=(lEESf"+ES^+EE 


dUi 


dUji 


dU. 


jk 




(29) 


Since 


9Uik _ dUik dU, 

/ J ( J f^p ■ ■ ' 


k^i 
we have 






drii ’ 


(30) 


Using this property, the total force on particle i can be 
derived as 



9Uj 

^ dri 


_ I ^Uj ^'f'ik 

^ I dri 

\k^j ■> 

/ dUi dUj \ 

^\drij dr^i)' 



dUi 

dri 


(37) 


which is equivalent to Eq. m, and the pairwise force is 
simplified to be 


pTersoff _ f ^Ui _ ^Uj \ 

~\dri, dr.J- 


(38) 




2 ^ drij 


Uij + Uji + ^Uik + Ujk) 1 . (31) 

k^i,j 


Erom this, a pairwise force between two particles can also 
be defined for the many-body Tersoff potential: 


^Tersoff _ 


1 d 


2 dri 


Uij Uji + / ^ {Uik + Ujk) 


k^i,3 


(32) 

The total force can be expressed as a sum of the pairwise 
forces 


( 33 ) 


One can check that Eq. (|38)) reduces to Eq. dni) in the 
case of two-body interaction. We also point out that our 
force expressions for the Tersoff potential are only seem¬ 
ingly different from other alternatives. There should be 
no ambiguity for the calculation of the total force on a 
given particle. However, different formulations may lead 
to different computer implementations. A crucial advan¬ 
tage of our formulation is that the total forces for indi¬ 
vidual particles can be calculated independently, which 
is desirable for massively parallel implementation. The 
numerical calculations presented in this work were per¬ 
formed by a molecular dynamics code implemented on 
GPUs using the thread-scheme in Ref. [^. However, a 
detailed presentation of the GPU-implementation of the 
Tersoff potential is beyond the scope this paper, which 
will be presented elsewhere. 
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Another advantage of our formulation is that the per- 
atom virial for the Tersoff potential takes the same form 
as for two-body potential: 


•CTT Tersoff 
VV 2 




-I Tersoff 


(39) 


which is unambiguously defined for periodic systemsi^ 
This might not be exactly equivalent to what has been 
implemented in LAMMPS, where the per-atom virial is 
calculated as 

-JEE ®®• (40) 

k^i,j 


Here, -F’gg represent the force compo¬ 

nents on particle i associated with the two-body part 
due to particle j, the three-body part due to particle j, 
and the three-body part due to particle fc, respectively. 
Although Eq. (l40l) and Eq. (l39l) may result in the same 
total virial tensor for the Tersoff potential, they may not 
be equivalent when used to compute the heat current and 
lattice thermal conductivity. We will present numerical 
results to compare them. 


E. Heat current for the Tersoff potential 


We now derive the heat current expressions for the 
Tersoff potential, using the potential decomposition given 
by Eq. ((35l) . Using Eq. (l37)) . the first term on the RHS 
of Eq. ([HI) can be written as 


i i j^i 



drji) 


Vi. (41) 


Using Eq. (1551) . the second term on the RHS of Eq. (IH|) 
can be written as 



i 


EE 

I 


dU, 


■ 


(42) 


From these two expressions, we get the following formula 
for the potential part of the heat current for the Tersoff 
potential: 


rTersoff 

^pot 



dvj. 



(43) 


Again, one can get rid of the absolute positions by 
rewriting the above formula as: 


yTersoff 

^pot 



dvj. 



(44) 


A less symmetric form can also be readily obtained: 


rTersoff 

^pot 



or equivalently. 


rTersoff 

^pot 



dvji 



(45) 


(46) 


Therefore, the potential part of the heat current for 
the Tersoff potential is not equivalent to the stress-based 
formula given by Eq. (155)) . One can check that, in the 
case of two-body interactions, the heat current expres¬ 
sions in Eqs. (I43II46I) for the Tersoff potential reduce to 
those for the two-body potential in Eqs. (HMni). 

Apart from the velocities Vi and relative positions r^, 
the only nontrivial terms in the force and heat current 
expressions are and the latter being able to 

be obtained from the former by an exchange of i and 
j. An explicit expression for the former is presented in 
appendix [^ 

In appendix |B1 we show that Eq. (H51) is equivalent to 
the one derived by Hardj^ at the quantum level for gen¬ 
eral many-body interactions. In the following, we refer to 
Eq. (H51) as the Hardy formula and Eq. (1501) as the stress 
formula. 

There has been some confusion about the seemingly 
different heat current expressions for the Tersoff potential 
in the literature. Guajardo-Cuellar et alt^ and Khadem 
et ali^ compared several expressionsi^d^dJrJ^iS^ in the 
literature. From their results, it seems as if all of these 
expressions were inequivalent. In appendix m we show 
that many of them are equivalent to the Hardy formula. 


F. Generalization to other many-body potentials 

Besides the Tersoff potential, the Brenner potentialii 
and the Stillinger-Weber (SW) potentiali^ are also widely 
used in the study of covalently bonded systems. Here, we 
first show that the derivations for the Tersoff potential 
can be generalized to these potentials and then summa¬ 
rize our results for a general many-body potential. 

The generalization to the Brenner potential is straight¬ 
forward. The pair potential Uij for this takes the same 
form as that for the Tersoff potential [Eq. (1511) ]. The 
bond-order function bij, hence Uij, is only a function of 
the position difference vectors originating from particle 
i, although the explicit form of bij in the Brenner po¬ 
tential is more complicated. This is the only property 
we used to derive the pairwise force expression [Eq. (1321) ] 
for the Tersoff potential. Therefore, the same pairwise 
force expression also applies to the Brenner potential. 
Using the same potential partition as for the Tersoff po¬ 
tential, Ui = 4 I can arrive at a simplified 

pairwise force expression [Eq. (I5H1) ] and the Hardy for¬ 
mula [Eq. (H51) ] of heat current, as in the case of the 
Tersoff potential. 
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We next consider the SW potential. The total poten¬ 
tial energy consists of a two-body part and a three-body 
part, the latter being given ass^^ 


+ hjki + hkij), (47) 

i j>i k>j 


where 

hijk = A exp 


Tij Qj T iiz d 


cos 9ijk + - ) . (48) 


Here, A, 7 , and a are parameters and cosOijk is defined 
as in Eq. (1251) . Similar definitions apply to hjki and hkij. 
It is clear that hijk is symmetric in the last two indices: 
hijk = hikj. Using this property, we can reexpress the 
three-body part of the total potential as 


{hijk + hjki + hkij), (49) 

i j^i kiiij 


which can be further simplified as 


In fact, the pairwise force formula and the Hardy for¬ 
mula of heat current apply to any many-body potential, 
because the crucial property we have used in the above 
derivations, i.e., the pair potential Uij (or the site poten¬ 
tial Ui) is only a function of the set of vectors {'i'ij}j^i, 
is satisfied by any empirical potential: any other posi¬ 
tion difference vector can be expressed as the difference 
of two vectors in this set. In other words, the vectors 
{rij}j:^i form a complete set of independent arguments 
for any pair or site potential associated with particle i. 
We can summarize our formulations as follows. For a 
general classical many-body potential, 

U = '£uA{njh^i), (56) 

i 

there exists a pairwise force between two particles i and 
h 


r - r 


dUi dU, 


f)v ■ ■ ^ 

U! 


(57) 


= (50) 

i j^i k^i,j 

Without referring to any potential partition, but noticing 
that hijk is only a function of the position difference vec¬ 
tors originating from particle i, one can derive a pairwise 
force expression for the three-body part: 


a well defined virial tensor for periodic systems, 

W = -^EE^b®^b> (58) 

i jiii 

and a well defined potential part of the heat current for 
periodic systems. 


p(3) _ p(3) 


(51) 


1 


dhi^ 


km 


E E 


^hjkm 


2 dr., 

= 


■ (52) 

With a definition of the site potential, 

^ 2 '^ = ^ E E with C/(3) = ^ I/f \ (53) 

ji^i k=^i,j * 


the above pairwise force expression can be simplified to 


.,3, ac/7 suf' 

drij dr,i 


(54) 


This is formally the same as that for the Tersoff potential, 
the only difference being the form of the site potential. 
Adopting the above potential decomposition, and notic¬ 
ing that Ui is only a function of the position difference 
vectors originating from particle i, one can conhrm that 
the potential part of the heat current also takes the form 
of the Hardy formula: 




The existence of a pairwise force for classical many- 
body potentials, albeit not surprising according to the 
principles of classical mechanics, has not been widely rec¬ 
ognized in the community. Without an explicit expres¬ 
sion for the pairwise force, much effort has been devoted 
to constructing general expressions for the virial tensor in 
periodic systemsi^i^^ Our formulations are thus not only 
useful for thermal conductivity calculations based on the 
Green-Kubo formula, but can also find application in the 
study of properties related to the stress tensor. 


III. VALIDATING THE PAIRWISE FORCE 
EXPRESSION 

One of the interesting results from the previous section 
is that the interatomic forces for empirical many-body 
potentials are in fact totally pairwise. Here, we take the 
Tersoff potential an example and present numerical evi¬ 
dence for the correctness of the pairwise force by compar¬ 
ing the forces calculated using Eq. dSZD and those using 
the finite-difference formula, 

U{--- ,ri- Axe,.,,-- -)-U{--- ,ri + Axe^,, ---) 


(60) 















IV. APPLICATIONS ON THERMAL 
CONDUCTIVITY CALCULATIONS 
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FIG. 1. (Color online) (a) Total forces in the ^-direction 
on individual carbon atoms in a configuration generated by 
randomly shifting the positions of all the atoms from the per¬ 
fect graphene structnre by a small amount. The small solid 
dots and the larger open circles represent the results by using 
Eqs. daZl) and (feol). respectively, (b) Force differences between 
those obtained by using Eqs. (EH) and IMl) as shown in (a). 
Note that the testing system has more than 50 atoms, but 
only the data for the first 50 atoms are shown for clarity. 
Results for the y- and 2 -directions are similar. 


Here, Fi^ is the x-component of the total force on particle 
i when the system is in a specific configuration and AxBx 
is a small displacement vector of particle i along the x- 
direction from its original position . One can similarly 
consider the other directions, but the results for all the 
directions are quite similar and we only present the re¬ 
sults for a single direction for simplicity. We have checked 
that the forces calculated by the finite-difference method 
do not change over a wide range of Ax. 


The results for the comparison are shown in Fig.[TJ The 
system corresponds to a graphene sheet perturbed from 
the perfect honeycomb structure by randomly shifting 
the positions of all the atoms by a small amount. One can 
see that the forces on each particle calculated by the pair¬ 
wise force expression Eq. (EZl) and the finite-difference 
expression Eq. (l60l) are practically the same, with the 
relative errors being as small as about 10“®. This com¬ 
parison thus confirms the correctness of the pairwise force 
for the Tersoff potential unambiguously. 


With the force expression validated, we are now in a 
position to apply the heat current formulations to study 
lattice thermal conductivities of various kinds of mate¬ 
rial. To be specific, we present results obtained by using 
the Tersoff potential, which has been applied extensively 
in the study of thermal transport properties of silicon, 
diamond, graphene, and CNT. The Tersoff parameters 
used for diamond and silicon are taken from Ref. 0 
and those for graphene and CNT are the optimized ones 
obtained by Lindsay and Broidoj^l To be specific, we 
only consider isotopically pure and ^®Si in our simu¬ 
lations, although our method is not limited to this case. 
When calculating the thermal conductivity of graphene 
and CNT, one has to specify the effective thickness of 
the graphene sheet. We have chosen it to be 0.335 nm. 
We use cubic simulation cells for silicon and diamond 
and roughly square-shaped simulation cells for graphene. 
The time step of integration in the MD simulations is 
chosen to be 1 fs for most of the simulated systems, 
but for smaller carbon systems, we found that smaller 
time steps are desirable. The evolution time in the equi¬ 
libration stage (canonical ensemble, where temperature 
is controlled) of the MD simulation lasts one to several 
nanoseconds, depending on the simulations cell size. The 
heat current data are recorded every 10 steps in the pro¬ 
duction stage (microcanonical ensemble, where temper¬ 
ature is not controlled). We only consider systems with 
zero external pressure and the lattice constants for sili¬ 
con at 500 K and diamond at 300 K are determined to 
be 0.544 nm and 0.357 nm. Eor grahene and CNT at 300 
K, the average carbon-carbon distance is determined to 
be 0.144 nm. 


A. Silicon 

We start presenting our results by considering silicon. 
Eigs. Ha-e) show the RTCs [given by Eq. ([T])] for silicon 
at 500 K with different simulation cell sizes N. Eor a 
given N, there are large variations between the indepen¬ 
dent simulations associated with different sets of initial 
velocities in the MD simulations. Despite the variations, 
a well converged RTC can be obtained by averaging over 
sufficiently many independent simulations, along with es¬ 
timations of an average value and the corresponding error 
estimate for the converged thermal conductivity. In this 
work, we determine them in the following steps (for a 
given N)-. 

1. Determine (by visual inspection) a range of cor¬ 
relation time where the averaged RTC has 

converged well. 

2. Calculate the average values of the RTCs for the in¬ 
dependent simulations over the range of correlation 
time determined in the last step. 









(147±2)W/m 

. .p T 

-K 

.1... 

--"f i-< 



calculations. The Green-Kubo formula is, in principle, 
only meaningful for infinite systems, i.e., systems in the 
thermodynamic limit. However, in practice, one can only 
simulate systems with finite simulation cell sizes, with 
periodic boundary conditions applied along the directions 
which are thought to be infinite to alleviate the finite- 
size effects in those directions. One can then check if the 
results converge with increasing simulation cell size. 

Figure mf) presents the converged thermal conductiv¬ 
ities of silicon at 500 K obtained by using different sim¬ 
ulation cell sizes: N = 512, 1000, 1728, 2744, and 4096. 
It can be seen that they do not show a systematical de¬ 
creasing or increasing trend with increasing N. 

Due to the small finite-size effects, we can take the 
average values of thermal conductivity for different sim¬ 
ulation cell sizes as independent simulation results and 
obtain an average value and the corresponding error es¬ 
timate. In this way, we obtain the final result, (147 ± 2) 
W/m-K, which is in good agreement with that obtained 
by Howellj^ (155±4) W/m-K. Note that Howell used the 
direct method with the same Tersoff parameters. This 
comparison thus further confirmed the equivalence be¬ 
tween the direct method and the Green-Kubo method, 
as has been shown by Schelling et al^ for SW silicon. 


FIG. 2. (Color online) (a-e) Running thermal conductivities 
as a function of correlation time for silicon with different sim¬ 
ulation cell sizes at 500 K. The thinner (and lighter) and the 
thicker (and darker) lines represent the results of independent 
simulations with different initial velocities and the ensemble 
average over the independent simulations, respectively, (f) 
Thermal conductivity as a function of the simulation cell size 
N. Markers with error bars represent the average values and 
the corresponding standard errors for a given N. The solid 
line indicates the average (147 W/m-K) over the 5 simula¬ 
tion cell sizes and the dashed lines indicate the corresponding 
standard error (±2 W/m-K). 


3. Take the mean value and standard error (standard 
deviation divided by a/M, where M is the number 
of independent simulations) of the average values 
obtained in the last step as the average value and 
error estimate, which are represented by an open 
circle and the corresponding error bar in Fig. [21(f) 
for a given N. 

To determine [ti, ^ 2 ], we have to ensure that the aver¬ 
aged RTG is sufficiently smooth. The smoothness can be 
enhanced by increasing either the simulation time of 
the individual simulations or the number of independent 
simulations Ng- More precisely, it is determined by the 
product Ngts- We found that a value of Ngts = 200 ns 
is enough for silicon at 500 K. It can be seen that all 
the averaged RTCs in Figs. |2ja-e) are rather smooth and 
[^ 1 ,^ 2 ] = [400 ps, 500 ps] is a fairly good choice for the 
converged time interval. 

Before comparing our results with previous ones, we 
need to further check possible finite-size effects in the 


B. Diamond 





FIG. 3. (Golor online) Same as Fig. O but for diamond at 
300 K. 


We next consider diamond. The RTCs at 300 K with 
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5 simulation cell sizes, N = 512, 1000, 1728, 2744, and 
4096, are shown in Figs. Ha-e) and the corresponding 
converged values are presented in Fig.[3l[f). The averaged 
RTCs converge earlier than those for silicon. Here, it can 
be seen that the converged time interval can be chosen 
to be [ti,t 2 ] = [150 ps, 200 ps]. Due to the shorter corre¬ 
lation time required for converging, the total simulation 
time required for obtaining smooth curves of the RTC is 
shorter than that for silicon, being about Ngts = 100 ns. 

As in the case of silicon, there is no systematical de¬ 
creasing or increasing trend with increasing N. Our cal¬ 
culated thermal conductivity averaged over the 5 simula¬ 
tion cell sizes is (1950 ± 40) W/m-K. Using the Brenner 
potentialii and the Green-Kubo method, Che et ah^ 
obtained a converged value of about 1200 W/m-K for 
isotopically pure diamond, which is about one third 
smaller than ours. This difference can be understood 
by noticing that the original Brenner potential is more 
anharmonic than the original Tersoff potential, as has 
also been noticed in the study of CNT and graphenei^l 
Experimentally, the thermal conductivity of isotopically 
pure diamond at room temperature is about 3000 
W/in-K^ larger than both of our results. The differ¬ 
ence between theoretical and experimental results may 
result from an excessive anharmonicity of the empirical 
potentials. 
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FIG. 4. (Color online) Same as Fig. O but for graphene at 
300 K. 


C. Graphene 

The above results are for 3D bulk materials. We now 
turn to study low-dimensional materials, first considering 
2D graphene. The RTCs at 300 K with 5 simulation cell 
sizes, N = 960, 3840, 8640, 15360, and 24000, are shown 
in Figs. UDa-e), with the corresponding converged values 
presented in Fig. HDf). For each N, a total simulation 
time of Nstg = 500 ns is required to obtain an average 
RTC well converged in the time interval of [U, f 2 ] = [250 
ps, 500 ps]. 

As in the case of diamond and silicon, the thermal con¬ 
ductivity of graphene does not increase with increasing 
simulation cell size. In fact, the contrary is true when N 
is smaller than 10“^, as found by Pereira and Donadio<^ 
Similar results have also been obtained by Zhang et ah^ 
for smaller N. The increasing of the simulation cell 
size has two opposite effects: (1) It allows more long- 
wavelength phonons, which can increase the thermal con¬ 
ductivity; (2) It also allows more phonon scattering, as 
suggested^^ by Ladd et aL, which can decrease the ther¬ 
mal conductivity. In 2D graphene, more phonon scatter¬ 
ing can be induced by the acoustic flexural modes with 
increasing out-of-plane deformation, which is positively 
correlated to the simulation cell sizei^ When the simu¬ 
lation cell size is relatively small, the second effect may 
dominate, resulting in a decreasing thermal conductivity 
with increasing simulation cell size. When the simula¬ 
tion cell size is relatively large, these two effects largely 
compensate each other, resulting in converged thermal 


conductivity with increasing simulation cell size. 

The thermal conductivity of graphene at 300 K av¬ 
eraged over the 5 simulation cell sizes is (2700 ± 80) 
W/m-K. Using the optimized Brenner potential^i and 
the Green-Kubo method, Zhang et reported a 

converged value of (2900 ± 93) W/m-K for graphene 
at 300 K, which is slightly larger than ours. This dif¬ 
ference may be explained by the fact that they have 
used smaller simulation cell sizes, which, according to 
the discussion above, results in larger thermal conduc¬ 
tivity for graphene. On the other hand, Haskins et al^ 
reported a value of 2600 W/m-K based on the Einstein 
formulation^ which is in good agreement with ours. 

It is interesting to point out that our estimate of the 
thermal conductivity for graphene at room temperature 
is compatible with NEMD calculations (ming the same 
Tersoff potential parameters) in Ref. [^, which give 
K ~ 2300 W/m-K with a simulation length of about 1.5 
fim. If we take the consistency between the Green-Kubo 
method and the NEMD method as granted, this compar¬ 
ison indicates that the NEMD results have not been con¬ 
verged up to a simulation length of 1.5 ^m. In fact, both 
the NEMD results and the experimental datai^^ suggest 
a logarithmic length-dependence of thermal conductiv¬ 
ity of graphene at the micrometer scale. On the other 
hand, whether the thermal conductivity is upper-limited 
or not in the infinite-size limit has been largely debated 
recently<^“— Our results provide evidence that the ther¬ 
mal conductivity of an extended (macroscopic) graphene 
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sheet is finite, although at the micrometer scale k still 
depends on the length of the graphene patch. 

D. (10, 0)-carbon nanotube 


Ref. [ 4 ^ employed the original parameter set provided 
by Tersofffi^; Ref. [i^ used the stress formula as imple¬ 
mented in LAMMPS, which also results in smaller values 
of K comparing with the Hardy formula, as we show be¬ 
low. 
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FIG. 5. (Color online) Same as Fig. O but for (10, 0)-CNT 
at 300 K. 


Last, we examine the longitudinal thermal conductiv¬ 
ity of CNT. To be specific, we consider a (10, 0)-CNT, 
without a detailed study of the effects of chirality and 
radius. The RTCs at 300 K with 5 simulation cell sizes, 
N = 2000, 4000, 6000, 8000, and 10000, are shown in 
Figs. EDa-e), with the corresponding converged values 
presented in Fig. EDf). For each N, a total simulation 
time of Ngts = 1000 ns is required to obtain an average 
RTC almost converged in the time interval of [^ 1 ,^ 2 ] = 
[500 ps, 1000 ps]. 

Compared with 2D graphene, the (10, 0)-CNT has 
even larger thermal conductivity: (3100 ± 68) W/m-K. 
This high value of thermal conductivity is mostly due 
to the long phonon wavelength (large phonon relaxation 
time) in QID CNTs^ as indicated by the slow conver¬ 
gence of K with respect to t. While there were debates 
on the size convergence of k for CNTs^“— our results 
do not suggest a divergent k with respect to the simu¬ 
lation cell length. Previously, the thermal conductivity 
for (10, 0)-CNT was calculated to be (1750 ±230) W/m- 
K in Ref. ^ (see also Ref. HI]) and (1700 ± 200) 
W/m-K in Ref. which are both smaller than the 

value obtained in this work, but due to different reasons: 


E. Comparing the stress and the Hardy formula 





correlation time (ps) 


FIG. 6. (Color online) Running thermal conductivities K,{t) 
as a function of correlation time for (a) silicon at 500 K, (b) 
diamond at 300 K, (c) graphene at 300 K, and (d) (10, 0)- 
CNT at 300 K obtained by using the Hardy formula (solid 
lines), the stress formula (dashed lines) and LAMMPS (dot- 
dashed lines). For each material, the line and the shaded area 
represent the averaged /t(t) and the standard error calculated 
from an ensemble of 10 independent simulations. 


Previously, we remarked that the stress formula [Eq. 
O] and the Hardy formula [Eq. (1461) ] are inequivalent 
for the Tersoff potential. Also, the per-atom virial as 
implemented in LAMMPS [Eq. (HOI) ] is not likely to be 
equivalent to ours [Eq. (l39l) ]. which would result in dif¬ 
ferent heat currents based on the stress formula. Here, 
we show the nonequivalence numerically. 

Eigure El shows the RTCs of (a) silicon at 500 K, (b) 
diamond at 300 K, (c) graphene at 300 K, and (d) (10, 0)- 
CNT at 300 K calculated using the Hardy formula, the 
stress formula in our formulation, and the stress formula 
as implemented in LAMMPS. We note the following ob¬ 
servations based on Pig. El 

(1) Por 3D diamond and silicon, all the three methods 
result in comparable results. 

(2) For 2D graphene, the RTC in the converged regime 
([250 ps, 500 ps]) obtained by the stress formula is 
about 1/2 of that by the Hardy formula and that by 
the LAMMPS implementation is about 1/3 of that by 
the Hardy formula. The LAMMPS results are consistent 
with previous oneSi ^^d^ 
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(3) For QID CNT, while the RTC in the converged 
regime ([400 ps, 600 ps]) obtained by the stress formula 
is comparable to that by the Hardy formula, that by the 
LAMMPS implementation is about 1/2 of that by the 
Hardy formula. The LAMMPS results are also consistent 
with previous ones 4^ 

From these observations, we conclude that the stress 
formula is generally inequivalent to the Hardy formula, 
and the LAMMPS implementation of the stress formula 
is also inequivalent to our implementation based on the 
pairwise force. Although we are not clear about the rea¬ 
son why the differences between these formulations are 
more significant in low-dimensional materials (especially 
2D graphene) than in 3D materials, our results can ex¬ 
plain an extraordinary low value of thermal conductiv¬ 
ity of graphene at 300 K, (280 ± 15) W/m-K, obtained 
by Mortazavi et al^ using LAMMPS and the (second- 
generation) Brenner potential.— Apart from the higher 
anharmonicity of this empirical potential compared with 
the optimized Tersoff potential, this small thermal con¬ 
ductivity could be attributed to the use of the stress for¬ 
mula implemented in LAMMPS. 


V. CONCLUSIONS 


In summary, we formulated force, stress, and heat cur¬ 
rent expressions of many-body potentials in MD simu¬ 
lations. After deriving these expressions for the Ter¬ 
soff potential in detail and briefly discussing their gen¬ 
eralizations to the Brenner potential and the Stillinger- 
Weber potential, we reached a set of universal expressions 
[Eqs. (I57II59P ] which apply to general many-body poten¬ 
tials. 

The pairwise force expression [Eq. (1571) ]. whose exis¬ 
tence is guaranteed by the principles of classical mechan¬ 
ics, has not been widely recognized in the community 
so far. We presented explicit numerical evidence of its 
correctness by comparing with finite-difference calcula¬ 
tions and demonstrated its importance in the construc¬ 
tion of a well defined virial tensor [Eq. (1551) ]. With a 
reasonable potential partition, we arrived at the Hardy 
formula [Eq. ([5911 ] for the potential part of microscopic 
heat current used in lattice thermal conductivity calcu¬ 
lations based on the Green-Kubo formula. Many of the 
seemingly different formulations of the heat current in 
the literature were demonstrated to be equivalent to the 
Hardy formula. 

We have implemented the formulations for the Ter¬ 
soff potential on GPUs and obtained orders of magni¬ 
tude speedup compared to our serial GPU implementa¬ 
tion. While the details of the GPU-implementation is 
beyond the scope of this paper, we have applied it to cal¬ 
culate systematically the lattice thermal conductivities 
of various kinds of material, including 3D silicon and di¬ 
amond, 2D graphene, and QID GNT, with emphasis on 
the effects of the simulation time and simulation cell size. 
We demonstrated the correctness of our formulations by 


comparing our results with previous ones. Last, we pro¬ 
vided explicit evidence to the nonequivalence between 
the Hardy formula and the stress formula as well as to 
the nonequivalence between the LAMMPS implementa¬ 
tion of the stress formula and our implementation based 
on the pairwise force. Particularly, we showed that the 
stress-based formulation underestimates the thermal con¬ 
ductivity of systems described by many-body potentials, 
and that this effect is more noticeable for low-dimensional 
systems. Our findings are very relevant for scientists 
modelling thermal transport in low-dimensional systems 
via molecular dynamics simulations. 
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Appendix A: Explicit expression for 


In this appendix, we present an explicit expression for 
which can be easily implemented in a computer 
language. 

Using the partition given by Eq. (ESD, we have 


du, _ idu,j , 1 ^ du,k 


dr 


I] 


2 dr 






dr,j ■ 


(Al) 


After some algebra, we have 
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where 


dr,. 


(A3) 


Since ^ = 0, we have 


m 

dr, ' 

I j^t 


jpot — X! X! 


(B7) 


d cos %fc 
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' ^3 L 


ik '^ij /I 

-^ cos 6 ^- 

r- • 

I ij 
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ijk 


(A4) 


and we have used the following notations: 

= dfA{r,j)/drij, fnirij) = dfR{r,j)/drij, 
fh^ij) = dfc{rij)/dr,j, b',^ = db,j/dCij, and 

9zjk = dg,jk/d cos 0ijk. 


Appendix B: Unifying different heat current 
expressions in the literature 

The derivation of the heat current expressions for a 
general lattice has been considered very early by Hardy 2 ^ 
at the quantum level. The potential part of the heat 
current was derived to be 


rHardy 

pot 




i j^i 


1 

' ih 


Pi 

2mi 




+ h.c., (Bl) 


where h is reduced Planck constant, and rrii are the 
momentum operator and mass for particle i, and h.c. 
stands for Hermitian conjugate. Using the identity 


[p„Uj] = 


(B2) 


the classical analog of Eq. m can be derived to be 
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Hardy 

pot 


= EE^. 


dr. 
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Using Eq. (IMl) . we have 

dUj drjk dUj 
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dr,}, dri dr,} ’ 


and 
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\ dr 

- / \ 72 
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(B3) 
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(B5) 


This equation is identical to Eq. (l46ll and thus equivalent 
to all the expressions in Eqs. (I43M5D . 

We now show that many of the seemingly inequivalent 
expressions of the potential part of the heat current for 
the Tersoff/Brenner potential are equivalent to the Hardy 
formula. 

We first consider the one used by Li et al.^ which 
takes the following form: 


jLi 

pot 


dEi 




(B 6 ) 


which has the same form as that used by Dong et al.r^ 
By noticing that [where we have used Eq. (IMl) ] 


dUi dUi dru, dU, 

dr ^ 




we have 
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i j^j^i 


dU, 
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(B 8 ) 


(B9) 


which is exactly Eq. (1451) and is thus equivalent to the 
Hardy formula. We also note that the one used by Berber 
et alr^ is exactly the Hardy formula. 

We next consider the one derived by Che et al.j^ which 
takes the following form: 




i 3 k I 


dru 


(BIO) 


Since 


dUki ^ ^ (Bii) 


dr 


we have 
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Che 

pot 




■ / drji 


= EE 


(B12) 


which is exactly the Hardy formula. 

The Hardy formula is also equivalent to a seemingly 
different one derived by Chen et al.,— which reads (The 
original expression in Ref. [13 contains a typo, which has 
been noticed by Guajardo-Cuellar et al.J^) 


yChen _ _ \ „ dUij 

pot — o / . / C b 


i j^i 


dr. 


-^EEE 


rik^^-Vk. (B13) 


2 drk 

By a change of indices {k o j), the second term on the 
RHS of the above equation can be written as 


1 dUik 

- oEE E 
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(B14) 


which, combining with the first term, gives [using 
Eq. (15^ ] 
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It takes the same form of Eq. (lB7l) and is thus equivalent 
to the Hardy formula. 


Recently, Guajardo-Cuellar et ah^ also derived an 
expression for the potential part of the heat current. 
They have used the equation rm ^ their 

derivation, which means that the force on particle i was 
taken to be Fi = This is only valid for two- 

body potentials, and as such it is not valid for the Tersoff 
potential. We thus do not expect that their expression is 
equivalent to the Hardy formula. 


parts: 
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^pot 


and 


-lEE 

% j/z 



+ E 


Vik 


djh 

dru 



TLi2 _ 

^pot “ 


z i/z 



dri 


Vi 


+ E '^ok 


dU,j 

dvk 



(B17) 

It can be shown that and ^ 

Jpot'^^/2 if one assumes the partition of potential energy 
given by Eq. (1551) . However, they have in fact chosen a 
different decomposition: 


= (B18) 


Last, we notice that Li et alwi^ also presented the po¬ 
tential part of the heat current as the sum of the following 


The calculated thermal conductivity is usually insensitive 
to the specific decomposition of the potential energy, as 
shown by Schelling et al^ for SW silicon. 
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